Waves and instability in a one-dimensional microfluidic array 
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Motion in a one-dimensional (ID) microfluidic array is simulated. Water droplets, dragged by 
flowing oil, are arranged in a single row, and due to their hydrodynamic interactions spacing between 
these droplets oscillates with a wave-like motion that is longitudinal or transverse. The simulation 
yields wave spectra that agree well with experiment. The wave-like motion has an instability which 
is confirmed to arise from nonlinearities in the interaction potential. The instability's growth is 
spatially localized. By selecting an appropriate correlation function, the interaction between the 
longitudinal and transverse waves is described. 
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I. INTRODUCTION 

Flows at a micron scale, called microfluidic flows, are 
of interest in fields such as molecular analysis, molecular 
biology and microelectronics [1]. Some microfluidic flows 
include a dispersed phase, such as droplets, bubbles, bio- 
logical cells or colloidal particles. In microfluidics, there 
is a flow in a channel that has at least one dimension that 
is hundreds of microns or smaller. In such a channel, a 
dispersed phase such as water droplets can be made to 
align in an array consisting of a single row, as sketched 
in Fig. 1; structures like these have been termed a 1- 
D array [2, 3], a 1-D stream [4], or a 1-D crystal [5-7]. 
Such 1-D arrays have applications in protein crystalliza- 
tion [2, 3] and flow cytometry [9], and they can exhibit 
physical phenomena such as waves [4-7] and instabilities 
[5-7, 10]. 

In this paper, we report a numerical simulation that 
mimics the experiment of Beatus et al. [5] . In that experi- 
ment, a 1-D array of water droplets was dragged by a flow 
of oil in a microfluidic channel, as sketched in Fig. 1. Due 
to the friction they experience on the floor and ceiling of 
the channel, the water droplets move more slowly than 
the oil. As they were injected into the oil flow, the wa- 
ter droplets at first self-organized in an equilibrium along 
the centerline of the channel, and they were positioned 
in a single row with a highly regular spacing, and the 
authors described this regularly-spaced array of droplets 
as a 1-D crystal. One discovery in the experiment of 
Beatus et al. was that the droplets fluctuate about their 
equilibrium positions with a wave-like motion and these 
oscillations are surprisingly not overdamped by viscous 
dissipation. This wave-like motion, termed phonon by 
the experimenters, included displacements in two direc- 
tions with respect to the flow, parallel and perpendicular, 
corresponding to longitudinal and transverse waves, re- 
spectively. Another finding in the experiment was an 
instability that developed farther downstream. In this 
instability, a local fluctuation in the 1-D array increased 
with time. As this instability grew, the displacements 
of droplets developed large amplitude fluctuation, which 
the experimenters attributed to nonlinear effects due to 
two factors: large amplitude motion and an interaction 



between longitudinal and transverse waves. Eventually 
the 1-D structure was destroyed by an extreme transverse 
displacement of droplets that allowed them to move past 
another in the longitudinal direction. The energy source 
for the instability was the flow of the oil through the 
channel [5]. 

Both the wave motion and instability in the dispersed 
phase were described by Beatus et al. [5-7] as the result 
of hydrodynamic interaction. One droplet disturbs the 
surrounding flow, which in turn disturbs another droplet 
via a drag force on the droplet due to the faster oil. Man- 
ifestations of this kind of hydrodynamic interaction have 
been reported for other microfluidic flows; these mani- 
festations include self-assembly patterns [11-13], shock 
waves [4, 7], and oscillatory waves. The latter were dis- 
covered experimentally by Beatus et al. [5, 7] and studied 
in simulations where the flow field scattered by droplets 
or particles is modeled as a superposition of the far-fields 
of dipoles [5, 7] or is calculated by a Stokesian-dynamics 
method [8]. 

Our simulation treats individual droplets as particles 
that move according to a well-defined equation of mo- 
tion. Instead of computing the entire flow field for the 
oil, we model the hydrodynamic interactions due to the 
flow field in the droplets' equation of motion. The model 
we use for the hydrodynamic interactions was developed 
theoretically by Beatus et al. [5, 7] who treated the oil 
flow as a potential flow. 

The simulation results we report here show agreement 
with the wave spectrum, dispersion relation, and insta- 
bility in the experiment of Beatus et al. [5]. We use the 
simulation to verify that the instability is due to non- 
linearity, as proposed by Beatus et al. [5]. We use the 
simulation results to describe the interaction of the lon- 
gitudinal and transverse waves by selecting an appropri- 
ate correlation function based on the time series of the 
longitudinal and transverse microscopic currents. 

The physics that will be studied here, in a 1-D mi- 
crofluidic array, is also relevant to other physical systems 
that have a single row of particles or atoms. These sys- 
tems include a line of colloidal microspheres confined by 
holographic optical traps [14], a line of charged micron- 
sized particles levitated in dusty plasmas [15, 16], and a 
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line of helium atoms adsorbed on bundles of single- walled 
carbon nanotubes [17]. In these examples, the particles 
require some mechanism to confine them so that they 
form a 1-D array, and the collection of particles can sus- 
tain wave-like motion, which is due to a combination of 
the confinement and the interaction between the parti- 
cles. 



II. SIMULATION 

We simulate a 1-D array of water droplets with the 
same configuration and parameters as in the experiment 
of Beatus et al. [5]. As sketched in Fig. 1, the channel is 
aligned in the x direction and it has a rectangular cross- 
section of width 250 /im and a height 10 /im that is the 
same as the droplet's radius R = 10 fim. The simulation 
begins after water droplets have been injected, so that 
they are initially spaced by a = 27 /im. In calculating 
the forces acting on the droplets, we assume an oil flow 
velocity of u™ u = 1730 /im/s, while the droplets move 
more slowly, at u,i = 380 /im/s. 

We numerically integrate the equations of motion of all 
the simulated droplets. A droplet is treated individually 
as a point object that responds to hydrodynamic forces, 
which are calculated to take into account of the finite 
radius R of the flattened droplet. The equation of motion 
for the n th droplet is [5] 

at U oil 

This equation of motion expresses a balance of two forces: 
the friction on the channel walls and the drag in the oil 
flow. The equation of motion in Eq. (1) includes only hy- 
drodynamic effects; other effects such as Brownian mo- 
tion due to the discreteness of molecules in the oil or 
water are neglected. The inertial term in Newton's sec- 
ond law is also neglected, so that unlike some equations of 
motion, Eq. (1) is a first-order differential equation. For 
each droplet, we integrate the equation of motion, so that 
the position r n and velocity dv n /dt of all the droplets are 
advanced simultaneously in a sequence of time steps. 
As in Ref. [5], we model the potential <f> as 

0(r) = + £ Mr-r m ), (2) 

where the first term on the right is a potential due to oil 
flowing uniformly at u^ u and the second term is a super- 
position of the potentials due to all other droplets m. An 
underlying assumption in Eq. (2) is a pairwise interaction 
of droplets, meaning that the droplet phase is sufficiently 
dilute that three-body interactions can be neglected. For 
such a dilute phase, it is appropriate to model the far- 
field interaction using only a dipole potential. As in [5], 
for the interaction potential <f>d we use 

Mr)=R 2 (uZ-u d )^, (3) 



which is the potential due to a dipole aligned with the 
flow in the x direction. Here, r = \/ x 1 + y 2 is the dis- 
tance from the droplet located at x = 0, y = 0. Equation 
(3) assumes that the droplet has been flattened by the 
channel to have the shape of a thin disc, and that the sur- 
rounding flow has a Poiseuille-like parabolic profile along 
the z axis and can be described in the xy plane as a 
potential flow, as in a Hele-Shaw cell. Equation (3) also 
assumes that the droplets are unconfined in the xy plane, 
i.e., we assume that the channel width is so great that its 
effects can be neglected. Such finite width effects could 
be accounted for by using a different potential [6, 7]. 

The force, which varies as — in the right-hand-side 
of Eq. (1), is nonlinear with respect to displacements of 
a droplet, when using Eqs. (2) and (3) to describe <fr. 
Most of our results are for this nonlinear case, although in 
Sec. Ill C 3 we will perform a test with a linearized form 
of — V(f>. The latter is obtained by performing a Taylor 
expansion of Eq. (3) for small displacements 5x and 5y 
from the equilibrium position x , where x = x + 5x and 
y = 5y. This expansion of Eq. (3) is 

x x \ x J 

- (*)" <« 

To turn off nonlinear effects in — V</>, we can retain only 
the first four terms in the brackets and neglect the higher- 
order terms (h.o.t.); we will do this only in Sec. IIIC3. 
For all other results, we will use the fully nonlinear force 
computed using Eq. (3) instead. 

Our simulation includes N = 256 droplets, and we 
use a periodic boundary condition in order to mimic an 
infinite 1-D array. The calculation is performed in the 
inertial frame of the droplets. To imitate the confine- 
ment provided by the floor and ceiling of the channel, the 
droplet is constrained to have no movement in the verti- 
cal direction, so that the equation of motion is integrated 
in only two coordinates, x and y. The simulation runs for 
a duration of 10 s. To advance a droplet's position, we 
use a fourth-order Runge-Kutta integrator with a fixed 
time step [19]. The time step was selected as 5.6 x 10~ 5 s 
by performing a test requiring that discrepancies in the 
kinetic energy of droplets are never significant over the 
simulation's duration. 

Because the droplets in the 1-D array are allowed to 
move in two directions, perpendicular and parallel to the 
flow, the 1-D array can sustain two modes, longitudinal 
and transverse, respectively. Due to the finite number 
of droplets, small amplitude motion can be decomposed 
into 256 discrete sinusoidal modes for the longitudinal 
motion, and the same for the transverse motion. Thus, 
the allowed values of ka in our simulation are ±7r/128, 
±27t/128, . . . , ±7r, where a = 27/im is the equilibrium 
spacing (lattice constant), and k = 2n/X is the wave 
number for a mode of wavelength A. 

Our initial conditions at t = are chosen to mimic 
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small random displacements of the droplets from their 
equilibrium positions. We start a droplet m at a posi- 
tion x rn = (m — N/2 + 5 x )a and y rn = 5 y a. Here, 5 X 
and 5 V are random numbers, with a mean of zero and a 
standard deviation of 0.01. The random numbers, which 
were different for each droplet, were chosen so that there 
would initially be equal spectral power in each of the 256 
longitudinal and 256 transverse modes. (This was done 
by initially giving each mode the same amplitude, but 
with a random phase.) 

As an approximation that is typical for molecular dy- 
namics simulations, we cut off the potential at large ra- 
dius, so that <f>d is given by Eq. (3) for r < r cut , or it 
is zero for r > r cut [18]. To select our cutoff radius 
r cut = 64a, we performed a test requiring that the ki- 
netic energy never develop a significant discrepancy. 

After running the simulation and obtaining time se- 
ries of the positions and velocities of all droplets, we use 
time scries delta for the longitudinal position x m and the 
velocity (v x , m 



T ) for each droplet to calculate 



and 



Xl(M) = V x,m{t)e 



j T (k,t) = %,m(*)e" 

m=l 



ikx m {t) 



ikxm (t) 



(5) 



(6) 



Equations (5) and (6) describe the fluctuating velocities 
of droplets in the longitudinal (x) and transverse (y) di- 
rections, respectively, for a specified value of the wave 
number k. In the literature for the physics of liquids, ji, 
and jt are sometimes called "microscopic currents" [20] , 
which can characterize the oscillatory wave-like motion 
of molecules in a liquid. Here, instead of molecules, we 
track the motion of droplets. We compute these currents 
for all 256 longitudinal and 256 transverse modes that 
are allowed for our simulation. 



III. RESULTS 

A. Droplet positions and velocities 

In Fig. 2 we present raw simulation data consisting of 
positions recorded at three times. At t = the simula- 
tion's initial conditions are shown, with droplets located 
almost at equal spacing a along the centerline of the chan- 
nel, with small random displacements in both x and y, 
as described in Sec. II. At t = 5 s, the displacements 
of the droplets have become large enough to be easily 
detected by our diagnostics, although they are still too 
small to be identified in a visual inspection of Fig. 2. At 
the end of the simulation, t = 10 s, the displacements 
of the droplets have grown more, as large as a, due to 
the instability. At i/a ~ —94 in Fig. 2(c), we observe a 
fluctuation resembling the 1+3 arrangement of droplets 



observed in the experiment [5, 7] and the simulation of 
Beatus et al. [7]. After 10 s, the droplets move past one 
another, i.e., change their sequence, and at that point in 
the development of the instability, the 1-D lattice struc- 
ture is essentially broken. The three panels shown in 
Fig. 2 are representative of the entire simulation movie, 
which we provide in the Supplemental Material [21]. 

We will use the raw data for droplet positions and 
velocities, and the currents computed from them using 
Eqs. (5) and (6), as the inputs to several diagnostics to 
characterize the waves and instability, as described next. 



B. Waves 

1. Wave spectra 

In Fig. 3 we present the wave spectra, which character- 
ize wave-like motion in the 1-D array of water droplets. 
These spectra were prepared by computing fast-Fourier 
transforms (FFT) of the currents jL(k,t) and jT(k,t) 
with respect to time for each allowed k, and then plot- 
ting the spectral power (the square of the modulus of the 
FFT) as a function of k and frequency to. This spectral 
power is physically related to the kinetic energy of the 
droplets. 

We find that our wave spectra in Fig. 3 agree with 
the experimental spectra [5]. The most obvious differ- 
ence is that our spectra do not include an artifact of un- 
wanted fabrication defects in the channel that deflected 
the droplet motion in the experiment. 

a. Dispersion relation. For both the longitudinal wave 
in Fig. 3(a) and the transverse wave in Fig. 3(b), the 
spectral power is significant only along a curve in the uj 
and k parameter space. This variation of the resonance 
frequency with wave number indicates the wave disper- 
sion relation. 

For comparison, we also show in Fig. 3 the theoretical 
dispersion relation derived in [5] from a linearized equa- 
tion of motion for small-amplitude motion with a dipole 
interaction 
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(7) 



We calculated the dispersion relation curves for 
this linear theory using the experimental values of 
the droplet spacing a and sound velocity C s — 
(2ty 2 R 2 u ( i/3o 2 u'^ 1 )(u^° 1 — u d ), with no free parameters. 
We find that the dispersion relation for the linear theory 
accurately matches our simulation results. 

The two waves propagate in opposite directions. In 
the inertial frame of the droplets, the transverse wave 
propagates in the +x (i.e., in the same direction as the 
oil flow), while the longitudinal wave propagates in the 
—x direction (opposite to the oil flow). These directions 
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of propagation can be identified in Fig. 3 by noting that 
for small \ka\, the phase velocities v$ = cj/k of trans- 
verse and longitudinal waves are positive and negative, 
respectively. These wave directions are the same as was 
observed in the experiment [5]. 

We find that the power in our spectra is not equally 
partitioned among the possible values of wave number 
k. The power becomes concentrated at the wave number 
where the frequency is maximum, ka ~ ±7r/2. We will 
comment upon this wave number later. 

b. Feature not predicted by linear theory. Although 
the dispersion relations observed in our simulation agree 
with the linear theory, we also observe a weaker feature 
that is not predicted by the linear theory of Eq. (7). This 
feature is seen near ka « ±0.8 and u & ±10 s _1 , where 
it is marked #3 in Fig. 3(b). Because these modes are 
not expected from a linear theory, there must be some 
nonlinear effects involved. 

Beatus et al. [5] suggested that nonlinearity arose in 
their experiment due to two factors: large amplitude mo- 
tion and an interaction of the longitudinal and transverse 
waves. Their suggestion motivates us to consider what 
type of wave-wave interaction it could be. 

One possible type of nonlinear interaction is the gen- 
eration of a third wave due by nonlinear mixing of two 
other waves obeying the criteria 

cji±cj 2 =^3! (8) 

and 

k 1 + k 2 = k 3 . (9) 

Such a process is called "parametric decay" in plasma 
physics [22] and "three- wave" mixing in optics [23] , and 
we will use the latter terminology. An argument in favor 
of this three-wave mixing hypothesis is that the u and 
k of the longitudinal waves marked #1 in Fig. 3(a) and 
the transverse waves marked #2 in Fig. 3(b), add up to 
match those of the feature marked #3 in Fig. 3(b), as 
in Eqs. (8) and (9). However, a counterargument is that 
feature #3 does not lie on an allowed dispersion relation 
for linear waves, unlike the best known cases of three- 
wave mixing in plasma physics [22] . 

2. Wave correlation 

Having found in Fig. 3 a possible indication of non- 
linear coupling of the longitudinal and transverse waves 
in the spectra, we now seek another diagnostic to in- 
dicate any interaction or synchronization of these two 
waves. Using the longitudinal and transverse currents, 
jL{kr,,t) and jT(k T ,t), respectively, we calculate a cor- 
relation function 

r ii, i. \ (jL(k L ,t)j T (kT,t + T)) 

C L T(kL,k T ,T) = — = (10) 

\j(\3L{k L M 2 )(\jT{k T ,t)\ 2 ) 



In Eq. (10), t is a time delay, and (• ■ ■) is an average over 
time t. The time interval for the averaging was selected 
to be as large as possible, given the 10 s duration of our 
simulation. The correlation function Clt is a complex 
number; we will report its modulus 

|C| = |C iT (*i. ) *r,0)|, (11) 

as a measure of the strength of the correlation between 
the two waves. For steady conditions, a value |C| = 1 
would indicate perfect correlation (although we should 
note that the conditions here are not steady, due to the 
growth of the instability). Since we are investigating 
whether a longitudinal wave and a transverse wave are 
correlated, we present |C| as a function of the longitu- 
dinal wave number kr, and transverse wave number kr, 
Fig. 4. 

The results in Fig. 4 indicate a significant correlation 
of the longitudinal and transverse waves. This is evident 
from the dark regions in Fig. 4, which indicate a high 
correlation \C\. 

To help interpret the correlation data in Fig. 4, we 
present Fig. 5 to indicate the conditions of highest cor- 
relation. We find that high correlations occur for waves 
that obey a matching condition for wave number, which 
is either 

k L = -kr (12) 

or 

fez,o = kra ± (tt — e), (13) 

where e«l. 

We also find that the waves that exhibit high correla- 
tion have a matching condition for frequency 

LOT, = LOT- (14) 

To demonstrate this matching condition, we combine 
the dependence of our correlation on wave number and 
the theoretical dispersion relations, Eq. (7), which are 
graphed on the edges of the main panel of Fig. 5. In ex- 
amining Fig. 5, one can start with the open circle on the 
longitudinal dispersion relation, and then follow dashed 
line HI to the left, where it is seen to have high corre- 
lations with two wave numbers for transverse waves, as 
indicated by vertical lines VI and V2. All three of these 
modes (indicated by open and solid circles) have frequen- 
cies that are identical. This observation leads us to the 
frequency matching condition, Eq. (14). 

While our correlation function indicates a significant 
correlation, it does not explain the mechanism for this 
correlation. This is, of course, a common limitation of 
using correlation functions. 

3. Wave fronts 

To visualize the spatiotemporal characteristics of the 
wave fronts, we present in Fig. 6 a space-time diagram, 
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which is a plot of droplet velocity as a function of posi- 
tion x and time t. We prepared this space-time diagram 
by combining the time series for position x(t) and one 
component of the velocity, for example v x (t) , to generate 
contours of constant v x in the parameter space x vs. t. 
A darker shade in Fig. 6 indicates a higher droplet speed. 
These space-time diagrams are useful for characterizing 
the spatiotemporal development of the droplet motion 
over the course of the entire simulation. 

Features that can be seen in the space-time diagram 
include wave fronts for the longitudinal and transverse 
waves. These wave fronts appear as sloped stripes. Be- 
cause the longitudinal waves propagate in the — x direc- 
tion, while the transverse waves propagate in the +x di- 
rection, these two kinds of wave fronts are sloped oppo- 
sitely in Fig. 6. 

The observed phase velocity of the waves is found 
by measuring the slope of the wave fronts in the space- 
time diagram. We find that vj, = — 100 /im/s for the 
longitudinal and v$ = fOO /im/s for the transverse waves. 

C. Instability 

We now present our results for the instability. We will 
find that the instability grows nonexponentially, with a 
distinctive spatial localization. We will also confirm that 
the instability requires: (1) nonlinearity in the forces and 
(2) a coupling between longitudinal and transverse mo- 
tions. 



disturbance is highly localized, with a width of 3a to 4a, 
i.e., it includes only three or four droplets. This width is 
consistent with the experimental observation of a fluctu- 
ation of 1 + 3 droplets [5, 7]. 

We note that the spatial width of « 4a of this distur- 
bance corresponds to a wave number of ka » tt/2. In the 
wave spectra, this same wave number of tt/2 was found 
to have a concentration of spectral power. Thus, we can 
suggest that the concentration of spectral power that is 
seen in the wave spectra is due to the spatial localization 
of the wave's growth. 



3. Requirement of nonlinearity 

We perform a test to assess the role of nonlinearities in 
the instability's growth. We do this by comparing results 
for the simulation run with different expressions for the 
potential: for the fully nonlinear run we use Eq. (3), while 
for the linear run we use Eq. (4) and retain only the first 
four terms on the right hand side of Eq. (4) . 

Results in Fig. 7 show that there is no growth of kinetic 
energy of the droplets when nonlinearities are turned off. 
Thus, we can conclude that the instability requires non- 
linearities in the potential. This conclusion confirms the 
suggestion of the experimenters [5, 7] that the instability 
is essentially the result of a nonlinearity. 



1. Temporal growth 

A feature that can be seen in the space-time diagram 
of Fig. 6 is a growth trend for the amplitude of the waves. 
As time progresses, the velocity amplitudes of both the 
longitudinal and transverse waves grow, as indicated by 
the shading in Fig. 6 appearing darker at larger times. 

To reveal the scaling of the growth with time, we can 
also examine a time series of the kinetic energy averaged 
over the entire length of the 1-D array, which we present 
in Fig. 7. We find that the kinetic energy grows with 
time, but unlike some hydrodynamic instabilities such 
as the Ray leigh- Taylor instability, this one has a growth 
that is not exponential with time. 

2. Spatial localization 

In addition to its temporal dependence, we can also 
characterize the spatial dependence of the instability's 
growth. This can be examined in the space-time dia- 
gram, Fig. 6. We find that the instability does not grow 
uniformly along the entire length of the 1-D array, but in- 
stead it develops in the form of isolated fluctuations. An 
example of this spatial concentration is seen prominently 
at t > 8 s, as marked * in Fig. 6(b). This large-amplitude 



4- Requirement of transverse motion 

We perform another test to determine whether motion 
in the longitudinal direction can grow due to an insta- 
bility in the absence of transverse motion. In this test, 
we constrain the droplets to move only along the x axis 
by loading the simulation with non-zero initial displace- 
ments only in the x direction. Otherwise, the simulation 
in this test was the same as for our other fully nonlinear 
ones. 

Results, shown with the open circles in Fig. 7, indicate 
no growth when the droplets are constrained to move 
only along the x axis. This finding demonstrates that 
the instability requires transverse motion. 

Considering that all low-amplitude fluctuations can be 
decomposed as a spectrum of longitudinal waves for mo- 
tion along the x axis and transverse waves for motion 
in the y direction, it is reasonable to conclude that the 
instability involves a coupling between longitudinal and 
transverse waves. Moreover, combining this conclusion 
with the previous one for nonlinearity, we have confirmed 
the suggestion of Beatus et al. that the instability in- 
volves an interaction of the longitudinal and transverse 
waves. While we have not determined the exact nature 
of this interaction, one possibility is the three-wave mix- 
ing mechanism described in Sec. IIIB 1. 
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IV. SUMMARY 

We performed a numerical simulation to study waves 
and an instability in a 1-D array of water droplets in 
a microfluidic channel. The droplets were modeled as 
point objects that interact with each other via a hydrody- 
namic potential. The oscillatory motion of the droplets, 
as characterized in the inertial frame of the droplets, ex- 
hibits two waves: a longitudinal wave that is backward 
with respect to the direction of oil flow, and a transverse 
wave that is forward. The simulation spectra for these 
waves agree well with the previous experiment that we 
simulate [5]. The droplet motion increases with time; 
this instability varies nonexponentially with time, and 



its growth is spatially localized. We performed two tests 
that confirm that the instability requires (f) nonlineari- 
ties in the hydrodynamic potential and (2) an interaction 
between longitudinal and transverse waves. A possible 
candidate for this interaction is a nonlinear three-wave 
mixing, which is suggested by a feature seen in our wave 
spectrum. 
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FIG. 1: (Color online) Sketch of a 1-D array of water droplets dragged slowly by a faster flow of oil in a microfluidic channel, 
as in the experiment of [5] . 
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FIG. 2: Droplet positions at three times in the simulation, shown in the inertial frame of the array of droplets. Only a portion 
of a total of 256 droplets is shown here, (a) At the time t = 0, the initial positions of the droplets are specified. We then 
integrate the fully nonlinear equation of motion, Eqs. (l)-(3), to advance the droplet position and velocity each time step, (b) 
At t = 5 s, a droplet fluctuates about its equilibrium position, (c) At t = 10 s, the displacements of the droplets have grown 
as large as a, due to the instability. 
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FIG. 3: (Color online) Power spectra of longitudinal (a) and transverse (b) waves, with the logarithm of spectral power shown 
as a function of frequency and wave number. Spectral power is mostly concentrated in narrow ranges, which are centered on 
the dispersion relations for a linear theory, Eq. (7), shown as dash lines. The feature near ka sa ±0.8 and w ~ ±10 s _1 , marked 
#3 in (b), is not predicted by the linear theory, and is interpreted here as an indication of nonlinear effects, possibly three-wave 
mixing of the waves marked #1 and #2. 




FIG. 4: (Color online) Correlation strength C|, as defined in Eq. (11), for longitudinal /cl and transverse kr waves. We find a 
strong correlation, for certain combinations of ki, and kr- The contour scale has a threshold of 0.2 to suppress noise. 
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FIG. 5: (Color online) Interpretation of correlation data in Fig. 4. Wave numbers that exhibit high correlations in Fig. 4 obey 
matching conditions for wave numbers shown with a dot-dash curve with k^a, = —kxa and dash curves with fcta = £;Ta±(7r — e) 
where £«1. We also find a matching condition for frequency, by comparing to the dispersion relations, Eq. (7), which are 
plotted along the edges of the correlation graph. As an example of these matching conditions, we consider a longitudinal 
wave k' as marked with an open circle. Following the horizontal line HI across the diagram, strong correlation occurs for the 
transverse waves with wave numbers — k' and k' + tt — S, which have the same frequency as the longitudinal wave. 
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FIG. 6: (Color online) Space-time diagrams for velocities (a) v x and (b) v y . To make the space-time diagram easier to examine, 
we present the data for only a portion of the array of droplets, 146 < m < 196. The sloped stripes in (a) and (b) indicate that 
longitudinal and transverse wave packets propagate in opposite directions. Note that the wave growth is spatially localized and 
that the high-amplitude wave packets at t > 8 s include only three or four droplets at the location marked *. 
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FIG. 7: (Color online) Instability growth, as measured by the time series of the mean squared velocity, averaged over the entire 
length of the 1-D array. The simulation was run three ways for this test. The solid curve is for the fully nonlinear case using 
Eq. (3) when integrating the equation of motion, as in Figs. (l)-(5). The dashed curve is for the linear case using Eq. (4). The 
open circles are results for the fully nonlinear simulation with droplets constrained to move only longitudinally, i.e., in the x 
direction. We find that for the nonlinear case there is rapid nonexponential growth, but for the other two cases there is no 
growth, indicating that the instability requires nonlinearity and transverse motion 



